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EFFECTIVE  VISCOSITY  IN  THE  SIMULATION  OF  SPATIALLY 
EVOLVING  SHEAR  FLOWS  WITH  MONOTONIC  FCT  MODELS 


1.  INTRODUCTION 

Recent  papers  [1-6]  have  reported  results  of  numerical  simulations  of  subsonic,  spatially 
evolving  two-  and  three-dimensional  planar  shear  layers  using  monotonic  FCT  models.  The 
numerical  model  solves  the  time-dependent,  compressible,  invisdd  conservation  equations  for 
mass,  momentum,  and  energy  in  three  dimensions  in  order  to  examine  the  evolution  of  large- 
scale  coherent  structures.  The  equations  are  solved  numericsJly  using  a  fourth-order  phase- 
accurate  FCT  algorithm,  directional  timestep-splitting  techniques  on  structured  grids,  and 
appropriate  inflow  and  outflow  boundary  conditions  [1,2].  This  approach  has  been  shown  to 
be  adequate  for  simulating  the  moderate  and  high-Reynolds- number  vorticity  dynamics  in  the 
transition  region  of  free  flows  and  reproduces  the  large-scale  features  of  the  flow  observed  in 
the  laboratory  experiments,  e.g.,  the  asymmetric  entrainment  [3],  the  distribution  of  merging 
locations  [4],  the  spreading  rate  of  the  mixing  layers  [5,6],  and  the  basic  three-dimensional 
features  of  the  vorticity  dynamics  [6,7].  In  this  paper,  we  address  some  of  the  numerical  issues 
of  resolution  which  are  important  for  validations  in  the  studies  of  physical  mechanisms  in  these 
simulations. 

A  large-eddy-simulation  approach  is  used  to  compute  the  evolution  of  the  transitional  flow 
dynamics.  The  nonlinear  FCT  high-frequency  filtering,  combined  with  the  conservative,  causal 
and  monotone  properties  of  the  algorithm,  are  expected  to  effectively  provide  a  minimal  subgrid 
model  by  maintaining  the  large  scale  structures  while  numerically  smoothing  the  scales  with 
wavelengths  smaller  than  a  few  computational  cells.  In  this  framework,  the  small  residual 
numerical  viscosity  of  the  algorithm,  combined  with  unresolved  small  scale  convection  at  high 
Reynolds  numbers,  mimics  the  behavior  of  physical  viscosity. 

The  objective  of  this  paper  is  to  study  the  effective  numerical  diffusion  of  the  algorithm 
in  the  low-Mach-number  simulation  of  &ee  mixing  layers,  and  its  dependence  on  gridding  and 
free-stream  velocity  ratio.  A  byproduct  of  this  work  is  to  assess  the  gridding  requirements  for 
an  FCT-based  shear-flow  model  including  physical  viscosity.  In  Sec.  2,  we  introduce  the  main 
steps  of  the  FCT  algorithm  as  used  in  the  model,  and  discuss  the  problem  of  measuring  the 
residual  numerical  diffusion  of  the  scheme.  In  Sec.  3.1,  we  review  the  calculation  of  the  laminar 
spread  of  a  free  mixing  layer  based  on  boundary  layer  theory.  The  theoretical  results  are  used 

Manuscript  approved  October  11,  1991. 


1 


as  reference  for  the  effective  measurement  of  the  global  numerical  diffusion  of  the  model  in 
Sec.  3.2.  The  final  conclusions  are  given  in  Sec.  4. 
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2.  THE  FCT  ALGORITHM 


We  restrict  our  discussion  to  the  one-dimensional,  explicit,  fourth-order  phase-accurate 
FCT  algorithm  [8]  used  in  the  shear-flow  applications.  The  numerical  model  solves  a  system  of 
generalized  continuity  equations  of  the  form 


df  djvf) 
dt  dr 


=  h, 


(1) 


where  /  typically  represents  mass,  momentum,  or  energy  densities,  h  is  an  appropriate  source 
term  dependent  on  the  flow  variables  and  their  spatial  derivatives,  v  is  the  fluid  velocity,  and 
r  is  a  spatial  variable.  The  algorithm  consists  of  a  two-step  predictor-corrector  scheme  which 
ensures  that  the  conserved  quantities  remain  monotonic  and  positive  when  so  required.  First, 
it  modifles  the  linear  properties  of  a  high-order  algorithm  by  adding  diffusion  during  convective 
transport  to  prevent  dispersive  ripples  from  arising  when  sharp  discontinuities  are  present.  The 
added  diffusion  is  then  removed  in  the  antidiffusion  phase  of  the  algorithm  in  such  a  way  that 
the  residual  numerical  diffusion  is  minimal  while  preserving  the  monotonidty  and  positivity  of 
the  scheme.  The  algorithmic  details  for  the  solution  of  Eq.  (1)  are  discussed  in  general  elsewhere 
[9].  In  order  to  introduce  the  basic  ideas  it  suffices  here  to  restrict  the  discussion  to  the  case  of 
the  mass  density  equation  (/  =  p,  h  =  0),  case  for  which  the  steps  are  briefly  discussed  below. 


In  the  first  stage  of  the  algorithm,  the  convection  phase  involves 

/' = 4”’  - 


(2a) 


followed  by  the  diffusion  phase 


/y  =  /'  +  -  /J-’.)).  (2‘) 

where  denotes  the  variable  /  at  grid  point  j  and  timestep  n,  j  +  ^  denotes  the  midpoint 
between  grid  points  j  and  j  -f- 1, 

6t  and  Sr  are  the  integration  timestep  and  the  grid  spacing,  respectively,  and  [8] 
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(4) 


1  .  1  2 
-  6 

In  turn,  the  correcting  antidifussion  stage  consists  of 

(5) 

where  the  raw  antidifFusion  fluxes  Me  deflned  by 

=  Mi+j(/i+i  “  fj)y  (6) 

and  $  denotes  the  corrected  antidifFusion  flux.  Monotonicity  is  preserved  by  the  flux  correction 
by  ensuring  that  the  difFusion  compensation  generates  no  new  maxima  or  minima  in  the  solution 
and  existing  extrema  are  not  accentuated  [8).  In  practice,  this  is  enforced  by  defining  the 
corrected  fluxes  as  follows. 


~  Afo*  0,M»n  ~  •^(/i  ~  »  (7) 

with  S  =  sign(fj+i  —  fj),  and  |5|  =  1. 

In  the  schemes  used  in  the  FCT  shear-flow  models,  the  antidifFusion  coefficients  Me 
deflned  in  terms  of  a  difFusion  parameter,  Z?c,  by 

=  (8) 

where  De  1-  For  De  =  I,  the  scheme  reduces  to  the  fourth-order  phoenical  FCT  scheme  [8]. 
A  special  case  arises  in  the  absence  of  convection,  when  =  0.  In  this  case,  in  the  vicinity 
of  a  sharp  discontinuity  the  antidifFusion  flux  survives  intact  through  the  flux  corrector, 
and  cancels  the  corresponding  difFusion  flux  exactly  if  Dc  =  1,  in  which  case  the  discontinuity 
is  preserved  after  the  FCT  cycle.  Standard  linear  stability  analysis  shows  that  for  Dc  <  1 
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the  scheme  remains  fourth-order-phase-accurate,  and  the  amplitude  accuracy  becomes  second- 
order,  with  the  second-order  error-term  for  the  amplitude  proportional  to  (1  -  Dc).  Thus, 
by  choosing  Dg  slightly  smaller  than  unity  (typically,  Dg  =  0.997  -  0.999  in  the  shear-flow 
applications)  a  small  nonvanishing  residual  numerical  diffusion  can  be  provided  by  the  scheme. 
A  local  estimate  of  the  residual  numerical  diffusion  introduced  by  Dg  in  the  limit  v=0  has 
been  communicated  to  the  authors  [10],  atnd  is  briefly  described  in  what  follows.  Combining  (2) 
and  (5)  and  comparing  with  the  second-order,  central-diflerenced  version  of  the  one-dimensional 
diffusion  equation. 


^-P^  =  0, 

dt  dT^  ’ 


(9) 


we  can  identify 


or  equivalently,  the  diffusion  coefficient  is 

Alternatively,  we  can  rewrite  this  expression  in  a  form  useful  in  calculations  performed 
with  a  fixed  Courant  number,  c  =  6l(lvl  -J-  a)peo*/^^» 

v  =  \(\-  D.)£r(W±i)£2i,  (12) 

0  c 

where  a  is  the  local  speed  of  sound.  The  inversely-proportional  dependence  of  P  on  c  reflects 
the  fact  that  a  larger  amount  of  residual  numerical  diffusion  results  from  taking  a  lower  value 
of  c  and  correspondingly  shorter  timesteps. 

Equations  (11)  and  (12)  give  local  estimates  for  the  numerical  diffusion  of  the  algorithm  in 
the  absence  of  convection.  Standard  one- dimensional  tests  show  that  sharp  discontinuities  are 
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maintained  by  the  FCT  algorithm  and  remain  confined  within  3-5  cells  [e.g.,  ref.  11].  This  in¬ 
dicates  that  the  actual  residual  numerical  diffusion  of  the  algorithm  is  nonlinear  and  somewhat 
dependent  on  the  solution  because  of  the  flux  correction.  As  a  consequence,  even  if  we  enforce 
constant  Sr  and  c  in  the  simulations  —  in  which  case  £q.  (12)  suggests  that  the  concept  of 
numerical  diffusion  can  be  physically  meaningful  —  the  effective  (global)  numerical  diffusion  in 
multidimensional  spatially  evolving  shear  flow  simulations  cannot  be  estimated  using  straight¬ 
forward  extrapolations  from  one-dimension2Ll  local  results.  Thus,  we  need  to  obtain  practical 
measures  of  the  effects  of  numerical  diffusion  when  non-linear  convection  and  compressibility 
are  present,  and  when  we  have  shear  discontinuities  rather  than  contact  discontinuities.  The 
present  approach  focuses  on  the  initial  laminar  spread  of  a  step-function  velocity  profile  due  to 
viscous  (numerical)  diffusion  in  the  limit  of  low  Mach  numbers.  In  this  case,  we  compare  the 
results  of  the  simulations  with  the  known  incompressible  solution. 
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3.  LAMINAR  SPREADING  OF  THE  MIXING  LAYER 
3.1  Boundary  Layer  Theory 

We  consider  the  laminar  spread  of  a  mixing  layer  initially  defined  by  a  step-function  profile 
for  the  streamwise  velocity  and  zero  transverse  velocity.  The  flow  configuration  is  indicated 
schematically  in  Fig.  la,  where  the  laminar  spread  of  the  velocity  profile  (as  given  by  g(T})  for 
X  >  0)  can  be  calculated  under  certain  flow  conditions  using  boundary  layer  theory  [12].  In  this 
regime,  the  flow  can  be  considered  virtually  incompressible,  and  the  problem  involves  a  thin 
shear  layer  with  v/u  «  1.  Under  these  conditions,  the  pressure  gradients  are  negligible  and 
the  steady-state  incompressible  equations  describing  the  problem  can  be  written 


du  dv 
dx  ^  dy 


=  0, 


(13) 


du  du  d^u 

u-r~  —  1/ - 

dx  dy  dy^ 


=  0, 


(14) 


where  u  is  the  kinematical  viscosity.  The  boundary  conditions  for  the  velocities  in  the  steady 
state  problem  are  specified  by 


u(0,y>0)  =  Ui,  «(0,y<0)=U2,  t»(0,y)  =  0,  (15a) 

«(x,y  = +oo)  =  Ui,  w(x,y  = -t-oo)  =  0,  (156) 

u(i,y  = -oo)  =  Ua,  w(i,y  = -oo)  =  0,  (15c) 

with  Ui>Ux. 

Defining  the  similarity  variable 

»7  =  y(ii//U2)"^  (16) 

the  solution  can  be  expressed  in  terms  of  a  function  G{r))  related  to  the  stream  function  $  and 
the  velocities  through 
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4-  =  iuxU2)iG(v), 


(17) 


u 


dGjy) 

drj 


=  Uigir,). 


The  function  G{ti)  satisfies  the  ordinary  differential  equation 


(18) 


dr]^ 


drj^ 


which  must  be  solved  with  the  boundary  conditions 


G(0)  =  0;  G\-oo)  =  \  =  U2lUi;  G'(+oo)  =  1. 


(19) 


(20) 


Equation  (14)  is  solved  numerically  using  finite  differences.  Typical  solutions  obtained  for  the 
laminar  spread  of  an  incompressible  mixing  layer  are  shown  in  Fig.  lb  for  A  =  2,5, 10,oo. 
These  solutions  are  used  as  reference  to  measure  the  effective  numerical  (viscous)  diffusion  of 
the  algorithm  when  simulating  thin  shear  flows. 

3.2  Simulated  Mixing  Layer 

In  order  to  talh  about  a  meaningful  effective  numerical  diffusion  when  Dc  ^  I  and  measure 
its  effects,  we  perform  simulations  for  fixed  Courant  numbers  in  appropriate  laminar  flow  cases 
using  uniform  grids  in  the  regions  of  interest.  We  restrict  our  time-dependent  simulations 
to  the  limit  of  low  Mach  numbers.  In  this  limit,  the  flow  is  virtually  incompressible  and  a 
comparison  can  be  set  with  the  results  in  section  3.1  after  attaining  the  steady-state  regime  using 
a  time-marching  approach.  We  seek  a  measure  of  the  effective  viscous  diffusion  by  comparing 
the  numerical  spread  of  an  initial  step-function  streamwise  velocity  profile  with  the  laminar 
spread  predicted  by  boundary  layer  theory.  The  choice  of  $  1  is  dictated  by  the  interest  in 
numerically  simulating  the  large  scale  features  of  transitional,  free  shear  flows  in  the  limit  of  large 
(but  finite)  Reynolds  numbers,  while  maintaining  the  accuracy  of  the  scheme  as  close  to  fourth 
order  as  possible.  In  the  present  study  we  have  specially  focused  on  the  case  Dc  =  0.999,  which 
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has  been  the  choice  of  most  of  the  previous  shear-flow  simulations  using  the  FCT  numerical 
model.  The  dependence  of  the  results  on  is  then  discussed  at  the  end  of  the  Section. 

Shear  flows  are  highly  unstable,  and  in  order  to  compare  the  development  of  the  mixing 
layer  with  theoretical  results,  we  need  to  isolate  the  laminar  (viscous)  growth  of  the  mixing 
layer  from  the  growth  due  to  the  Kelvin-Helmholtz  instability.  An  unavoidable  initial  mismatch 
due  to  discretization,  between  flow  and  boundary  conditions  near  the  inflow  introduces  small 
transverse  velocity  perturbations  which  excite  the  Kelvin-Helmholtz  instability  in  the  shear 
layer  near  the  inflow  boundary  [3].  This  initial  instability  is  subsequently  followed  by  vortex 
roD-up  and  global  self- sustained  instabilities,  in  which  new  vortex  roll-ups  are  triggered  in  the 
initial  shear  layer  by  pressure  disturbances  originating  in  the  fluid  accelerations  downstream 
[13].  To  ensure  that  the  spread  of  the  streamwise  velocity  profile  is  due  only  to  the  residual 
numerical  diffusion  of  the  algorithm,  we  have  chosen  to  force  the  transverse  velocity  to  maintain 
its  initial  and  inflow  value  (zero)  throughout  the  computational  domain  during  the  simulations. 

Since  the  FCT  shear  flow  model  is  nearly  inviscid,  the  model  is  expected  to  be  meaningful 
for  high- Reynolds- number  transitional  flow  regimes  -  in  which  the  large-scale  flow  features  are 
independent  of  Reynolds-number  {Re).  The  approach  used  to  obtain  measures  of  the  effective 
residual  numerical  viscosity  involves  approximations  which  become  valid  in  this  limit  of  high 
Re,  where  v/u  =  0(i)  =  0(1/Re),  and  ^  is  the  thickness  of  the  mixing  layer.  This  approach 
is  also  consistent  with  the  thin-shear-layer  approximations  used  to  derive  Eq.  (14)  from  the 
incompressible  Navier-Stokes  equations.  In  particular,  as  in  boundary  layer  theory,  the  viscous 
diffusion  term  proportional  to  d^ujdy^  is  the  dominant  one  in  this  limiting  flow- regime,  and  also 
responsible  for  the  observed  numerical  spread  of  the  simulated  mixing  layer  in  the  framework 
of  the  present  study,  where  the  streamwise  gradients  are  much  smaller. 

With  an  explicit  unsteady  model  we  can  not  practically  afford  to  deal  with  very  small 
Mach  numbers  because  the  small  timesteps  dictated  by  the  Courant  condition  for  numerical 
stability  determine  very  long  corresponding  runs  to  reach  steady  state.  A  compromise  choice 
in  this  work  in  order  to  have  practical  calculations  and  at  the  same  time  nearly  incompressible 
flow,  is  to  use  mean  free  Mach  numbers,  A4  ~  0.3-0. 4,  for  which  the  compressibility  effects  as 
measured  by  the  relative  mass-density  variations  Ap/p  ~  Ad^/2  [14],  are  at  most  of  the  order 
of  4-8%. 
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Figure  2  shows  the  flow  conflguration  for  the  simulated  mixing  layer  consisting  of  two 
coflowing  streams  of  air  entering  a  long  chamber,  where  the  separation  between  walls  is  chosen 
sufficiently  large  to  ensure  that  the  mixing  layer  growth  is  unaffected  by  their  presence.  Inflow 
and  outflow  boundary  conditions  are  imposed  where  appropriate  in  the  streamwise  direction 

(x) ,  and  reflecting  free-slip  wall  boundary  conditions  are  required  in  the  transverse  direction 

(y) - 


The  inflow  and  outflow  boundary  conditions  were  developed  and  tested  for  multidimen¬ 
sional  FCT  calculations  [1-3].  The  inflow  boundary  conditions  specify  the  density  and  velocity 
of  the  jet  and  impose  an  homogeneous  Neumann  condition  on  the  energy.  These  conditions 
are  implemented  at  the  inflow  by  specifying  the  guard-cell  values  for  the  mass  and  velocity 
densities: 


PGi  —  Pioflow*  (21®) 

=  tiinflow  =  Uiy),  {21b) 

voi  =  0,  (21c) 

and  imposing  a  zero-gradient  condition  on  the  pressure: 

PG^=Pl,  (21d) 


where  the  subscripts  Gi  and  1  refer  to  the  guard  cells  and  to  the  first  row  of  ceUs  inside 
the  computational  domain  (at  the  inflow),  respectively.  The  guard-cell  values  of  the  energy 
are  calculated  through  the  equation  of  state  as  a  function  of  the  mass  density,  momenta  and 
pressure.  By  imposing  a  floating  condition  on  the  pressure  at  the  inflow,  finite  (unsteady) 
cross-stream  pressure  differences  are  allowed  to  appear  in  the  initial  shear  layer  in  response  to 
acoustic  waves  generated  by  events  downstream. 

The  boundary  conditions  at  the  open  boundary  downstream  approximate  the  time- 
dependent  flow  equations  at  the  boundaries,  by  linearizing  the  inviscid  flow  equations  and 
reducing  them  to  advection  equations  in  the  outflow  direction  (i), 

dQ/dt  -I-  uxo^dQldz  =  0,  (22) 
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where  uioc  is  the  local  z-component  of  the  velocity  near  the  boundary.  This  equation  is  dis¬ 
cretized  by  a  first-order  upwind  scheme.  The  result  is  a  relation  between  the  outflow  guard-cell 
value  as  a  function  of  and  , 

+  (23) 

where  n  and  n  ~  1  denote  the  current  and  previous  integration  cycles,  and  N  corresponds  to 
the  boundary  cell  [2].  In  addition, 

€  =  tiiocAt/AiN,  (24) 

At  is  the  integration  timestep,  and  Az^v  is  the  size  of  the  computational  cell  at  the  outflow 
boundary.  The  outflow  conditions  specified  by  Eq.  23  are  imposed  on  the  mass-  and  momentum- 
density  variables.  Information  about  how  the  flow  relaxes  to  ambient  conditions  is  provided 
by  specifying  the  ambient  pressure  Pamb  ^'t  infinity  and  imposing  a  relaxation  rate  on  the 
pressure  [1,3].  The  guard-cell  value  of  the  pressure  is  defined  through  the  expression  Pcf,  = 
Ps  +  {Pamb  —  Ps)  X  ^xsIxgn'  This  defines  Pc„  by  interpolating  between  the  pressure  value 
at  the  N  -th  cell  and  the  value  Pamb  specified  at  infinity.  The  guard-ceU  values  of  the  energy 
are  then  calculated  through  the  equation  of  state  as  a  function  of  the  other  flow  quantities. 
These  inflow  conditions  allow  the  pressure  at  the  inflow  to  vary  in  response  to  acoustic  waves 
generated  by  events  downstream.  An  important  result  of  allowing  this  feedback  to  occur  is  that 
physical,  self-sustained  global  instabilities  can  occur  naturally  in  the  calculations  [13]. 

Since  the  numerical  diffusion  depends  on  the  grid  spacing,  the  latter  is  kept  constant  in 
the  region  of  shear  layer  development.  Twenty  or  more  evenly  spaced  computational  cells  are 
used  in  the  transverse  direction  in  the  neighborhood  of  the  center  of  the  shear  layer,  which  is 
initially  defined  by  a  (one-cell)  step-function  discontinuity.  Other  cases,  involving  nonsquare 
uniform  griddings  and  geometrical  stretching  in  the  streamwise  direction  are  also  considered 
in  order  to  specifically  investigate  the  effects  of  non-uniformities.  The  particular  features  of 
the  computational  grids  used  in  this  work  are  summarized  in  Fig.  3.  The  actual  dimensions 
of  each  domain  in  Fig.  3  are  chosen  in  such  a  way  to  ensure  computational  efficiency  while 
allowing  an  adequate  description  of  the  evolution  of  the  shear  layers.  The  timesteps  used  in 
the  calculations  are  determined  from  the  free-stream  conditions  by  imposing  Courant  numbers 
in  the  range  0.4-0.46. 
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We  examine  a  number  of  different  cases,  in  which  we  study  the  spreading  of  u{y)  from  an 
initial  step-function  velocity  profile  at  z  =  0  as  a  function  of  streamwise  distance  x,  in  terms  of 
the  free-stream  velocity  ratio,  gridding,  and  Courant  number.  In  each  case,  the  time-evolution 
of  the  calculations  is  pursued  long  enough  to  ensure  that  the  initial  transients  flow  out  of  the 
computational  domain  and  a  steady  state  regime  is  attained.  We  reduce  the  calculated  profiles 
of  normalized  streamwise  momentum  g  =  iP‘fi)/iPoU2)  at  different  streamwise  locations  with 
the  similarity  variable,  rj  =  y(xv/U2)~^ ,  where  p  »  po  =  ambient  mass  density,  and  v  is  an 
adjustable  viscosity  parameter.  The  effective  numerical  diffusion  Vg  is  defined  to  be  equal  to  the 
viscosity  parameter  giving  the  best  fit  of  p(q)  to  the  laminar  mixing  layer  solution  g{T))  =  u/U2- 
The  fit  is  based  on  the  least  rms  deviations  <7i  and  C2  defined  by 


aliv)  = 


1 

N.Ny 


51 2  [j( »?'.••) 

i=i  «=i 


2 


(25) 


<t|(v)  = 


(Nr  -  l){Ny  -  1) 


tel  tel 


where  rm  =  T]i,i{v)  =  (yi,j/2.— »yiv,]  and  [xi.X2,—M*yv;]  ««  the  fixed  stream- 

wise  and  cross-stream  sampling  locations,  respectively,  and  A  is  the  difference  operator, 
A^(t)  =  n{i  +  1)  —  H{i).  This  approach  is  intended  to  obtsun  a  quantitative  global  mea¬ 
sure  of  the  effective  numerical  diffusion  of  the  algorithm  and  to  examine  the  extent  to  which 
the  numerically  simulated  profiles  can  be  reduced  to  a  similarity  solution.  The  procedure  is 
expected  to  give  an  estimated  effective  Reynolds  number  associated  with  the  (small)  scales  of 
the  order  of  a  few  computational  cells. 


Typical  results  corresponding  to  simulations  with  A  =  5  in  grid  o  of  Fig.  3  are  shown  in 
Fig.  4.  Figure  4a  shows  the  behavior  of  <t\  and  <T2  as  a  function  of  v.  Associated  with  the  par¬ 
ticular  grid  chosen,  the  figure  shows  least  deviations  for  v  =  u*  w  0.16cm*sec“*,  corresponding 
to  air  viscosity  at  the  standard  temperature  and  pressure  conditions  of  the  simulations.  Profiles 
of  g(T))  are  compared  with  g{T})  in  Fig.  4b,  where  the  markers  correspond  to  actual  grid  points 
and  distinguish  between  different  streamwise  stations.  The  spreading  of  the  initial  step-function 
profile  is  essentially  confined  within  3-5  computational  cells  for  the  range  of  streamwise  loca- 
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tions  considered,  and  the  calculated  points  collapse  quite  well  in  the  vicinity  of  the  theoretical 
curve.  The  optimal  fit  in  this  case  leads  to  an  effective  Reynolds  number,  Re^tt  =  Ublv^  »  1800, 
where  U  -  (C^i  +  I/2)/2,  and  6  =  6y,  is  taken  to  be  the  initial  thickness  of  the  shear  layer.  Based 
on  Eq.  (12),  with  =  2.0  x  cm/sec,  o  =  3.5  x  10*  cm/see,  and  c  =  0.5,  the  local  estimate 
of  the  numerical  diffusion  V  for  this  case  turns  out  to  be  nearly  three  times  larger  than  Ug. 
A  larger  local  prediction  for  the  effective  numerical  diffusion  is  expected,  since  Eq.  (12)  will 
estimate  the  transverse  numerical  diffusion  of  the  shear  without  accounting  for  the  streamwise 
convection  of  the  fluid  elements. 

The  numerical  diffusion  is  slightly  different  on  each  side  of  the  shear  layer.  This  can  be 
observed  in  Figure  4b  for  A  =  5,  which  suggests  that  the  diffusivity  is  smaller  on  the  slower 
side  (q  >  0).  The  dependence  of  this  diffusivity  difference  as  a  function  of  free-stream  velocity 
ratio  A  can  be  obtained  by  examining  the  partial  rms  deviations  <71+  and  <7i_,  corresponding 
to  q  >  0  and  t?  <  0,  respectively,  defined  by 


<^i±(v)=jy~  ’  (27) 

*  * 

In  Fig.  5  we  plot  <7i,  <ri+,  <ri_  for  the  extreme  values  of  velocity  ratio  considered  here,  namely 
for  A  =  2  atnd  A  =  oo.  The  separation  between  the  minima  of  (7i+  and  <ti-  can  be  used 
to  measure  the  difference  between  the  numerical  diffusivities  on  each  side  of  the  shear  layer. 
These  differences  are  less  pronounced  for  the  case  A  =  2  with  more  nearly  similar  free-stream 
velocities,  and  can  be  attributed  to  a  small  streamwise-velocity  dependent  contribution  to  the 
global  numerical  diffusion,  associated  with  the  integration  in  the  z— direction.  The  origin  of  this 
contribution  can  be  understood  locally  by  noting  that  the  antidiffusion  will  never  cancel  the 
diffusion  exactly  at  the  end  of  the  FCT  cycle  if  convection  cannot  be  neglected  (cf.  Eqs.  (4)  and 
(8)).  A  local  analysis  similar  to  that  in  Sec.  3.1  shows  the  presence  of  an  additional  velocity- 
dependent  diffusion  term  with  expected  behavior  of  the  form  0{cSxu^),  which  is  consistent 
with  the  observed  (global)  <lifferences.  Alternatively,  velocity-dependent  discrepancies  between 
numerical  and  theoretical  results  can  also  be  expected  to  some  extent  due  to  the  fact  that  the 
flow  is  solenoidal  in  the  theoretical  (incompressible)  case,  whereas  the  regime  of  the  compu¬ 
tations  is  characterized  by  low,  but  finite,  Mach  number.  More  specifically,  in  the  framework 
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of  the  present  computations,  dvfdy  is  independent  of  dujdx  and  vanishes  by  construction,  so 
that  the  computed  velocity  field  has  a  nonzero  divergence,  V.u  =  dufdx  ^  0.  Thus,  zones 
of  positive  (negative)  divergence  appear  during  the  development  of  the  mixing  layer,  as  the 
top  (bottom)  stream  accelerates  (decelerates),  which  could  also  explain  the  larger  (smaller) 
diffusivity  observed  in  the  computed  profiles. 

Figure  6  shows  the  dependence  of  the  rms  deviations  <Ti  on  the  free-stream  velocity  ratio  A 
in  grid  a  of  Fig.  3  for  fixed  U2  =2x10*  cm  sec~^.  The  optimal  value  Ve  is  not  very  sensitive  to 
changes  in  A,  which  is  consistent  with  the  local  estimate  for  the  numerical  diffusion  (£q.  (12)), 
since  (|v|  +  a)pfak  is  the  same  for  aU  cases. 

Figure  7  shows  the  typical  linear  dependence  of  on  the  diffusion  parameter  De,  in 
qualitative  agreement  with  the  linear  dependence  on  (1  —  Dg)  predicted  by  the  local  estimate 
(£q.  (12)),  but  with  a  slope  significantly  smaller  (e.g.,  of  the  order  five  times  smaller  on  grid 
a,  for  A  =  5).  As  pointed  out  above,  a  smaller  actual  effective  diffusivity  than  that  predicted 
by  £q.  (12)  can  be  attributed  to  the  fact  that  the  latter  accounts  for  the  local  cross-stream 
transverse  numerical  diffusion  of  the  mixing  layer  when  no  streamwise  convection  is  present. 

Figure  8  shows  a  comparison  similar  to  that  in  Fig.  4b  (Dg  =  0.999),  for  the  case  Dg  = 
0.990,  corresponding  to  the  highest  value  of  v,  in  Fig.  7.  The  spreading  of  the  step-function 
is  now  significantly  faster,  i.e.,  involves  more  computational  cells  at  shorter  distances  from  the 
ori^n  than  those  conseidered  in  Fig.  4b.  The  results  also  indicate  that  an  improved  s^eement 
with  the  theoretical  spread  is  attained  as  the  effective  numerical  diffusion  Vg  increases. 

Figure  9  shows  the  effect  of  grid  spacing  on  the  rms  deviations  for  A  =  5,  by  considering 
a  finer  (Fig.  3b)  and  a  coarser  grid  (Fig.  3c)  than  that  used  for  the  cases  in  Fig.  5  (grid  a 
in  Fig.  3),  with  spacings  twice  as  large  and  one-and-a-half  times  smaller,  respectively.  The 
observed  growth  of  v*  in  Fig.  7  implies  a  nonlinear  dependence  on  grid  spacing  6x, 


V,  -  {6x)°\ 


(28) 


where  at  =  ai{Dg)  decreases  from  qj  =  1.48  to  Qi  =  1.16  when  Dg  is  reduced  from  0.999 
to  0.990,  approaching  the  value  qj  =  1  (predicted  by  £q.  (12))  as  v*  increases.  Thus,  in 
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calculations  performed  with  fixed  c,  and  when  the  initial  velocity  profile  is  defined  by  a  smooth 
function  in  such  a  way  that  the  initial  shear  layer  thickness  is  grid-independent,  Eq.  (28)  implies 
an  increase  in  effective  Reynolds  number  with  grid  refinement  {6x  — ►  0)  according  to 


iZCeff  ~ 


(29a) 


In  the  case  of  the  initial  step-function  velocity  profile  considered  in  this  work,  for  which 
the  initial  thickness  of  the  shear  layer  is  of  order  Sx,  we  get 


Re^ft  -  (296) 

whereas  £q.  (12)  implies  a  value  independent  of  grid  size  for  Rctg.  The  dependence  of  the 
effective  numerical  diffusion  on  the  Courant  number  is  shown  in  Fig.  10  for  grid  c  in  Fig.  3. 
The  results  now  suggest  the  nonlinear  dependence 


u.  ^  (30) 

where  oj  =  a2(i?c)  increases  from  aj  =  0.49  to  03  =  0.83  when  Dc  decreases  from  0.999  to 
0.990,  in  contrast  with  the  inversely  proportional  dependence  on  Courant  number  in  Eq.  (12). 
As  noted  above,  the  functional  behavior  predicted  by  Eq.  (12)  (03  =  1)  is  approached  for  larger 
values  of  v^. 

The  effects  of  non-uniformities  on  the  gridding  are  examined  in  Figs.  11-13  for  Dc  =  0.999. 
We  first  examine  the  effects  of  changing  the  aspect  ratio  6x/6y  of  the  grid  while  maintmning  it 
evenly  spaced  in  each  direction.  The  plot  of  rms  deviations  in  Fig.  11  suggests  that  changing 
from  6x/6y  =  1  (grid  a  of  Fig.  3)  to  6x/6y  =  2.5  (grid  d  of  Fig.  3)  does  not  significantly  affect 
the  nature  of  the  effective  numerical  diffusion.  This  is  physically  expected  to  some  extent. 
It  indicates  that  for  aspect  ratios  6x/6y  ^  1,  v*  is  essentially  determined  by  the  (numerical) 
viscous  term  proportional  to  d^u/dy^  introduced  by  the  algorithm.  In  Fig.  12  we  include 
results  obtained  on  grids  d  and  e  of  Fig.  3,  both  having  a  basic  gridding  aspect  ratio  2.5  which 
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is  maintained  throughout  the  region  of  interest  in  one  case,  and  geometrically  increased  in 
a  significant  portion  of  the  domain,  in  the  other.  As  can  be  expected,  the  results  indicate 
that  if  much  larger  non-uniformities  in  the  gridding  are  introduced  through  grid  stretching,  the 
numerical  term  associated  with  d^ufdx^  can  no  longer  be  considered  negligible  and  the  effective 
diffusion  can  consequently  become  significantly  larger  than  in  the  uniform  case. 

In  order  to  further  assess  the  dependence  of  the  effective  numerical  viscosity  on  the  gridding 
aspect-ratio,  we  examine  the  results  of  planar  shear  flow  simulations  also  initialized  with  a  step- 
function  velocity  profile,  but  for  which  the  zero-transverse-velocity  restriction  is  not  enforced, 
so  that  vortex  roll-up  can  now  take  place.  Figure  13  shows  typical  instantaneous  visualizations 
using  contours  of  vorticity  (Q)  and  of  a  passive  scalar  (<f>)  convected  with  the  flow  velocity, 
where  4>  is  defined  to  be  initially  zero  for  the  slower  stream  and  unity  for  the  faster  stream.  The 
calculations  are  performed  on  grids  /  and  g  of  Fig.  3,  and  the  flow  is  organized  by  adding  a  small 
time-dependent  perturbation  to  the  streamwise  velocity  at  the  inflow  with  frequency  equal  to 
that  of  the  most  unstable  mode.  The  figure  shows  the  presence  of  vortex  rolls  resulting  from 
the  evolution  of  the  nonlinear  Kelvin-Helmholtz  instability.  The  frames  shown  are  associated 
with  the  same  physical  time,  and  identical  conditions  except  for  the  uniform  gridding,  which 
involves  aspect  ratio  6x/6y  =  1  in  one  case  (Fig.  13a),  and  6x/6y  =  2.5  in  the  other  (Fig.  13b). 
In  both  simulations,  the  Courant  number  is  c  %  0.4, 6y  =  0.015  cm,  and  A  =:  5,  with  air  streams 
under  the  same  temperature  and  pressure  conditions  as  in  the  cases  discussed  previously.  In 
particular,  the  conditions  of  the  simulations  are  such  that  the  expected  numerical  dlffusivities 
in  the  case  with  SxfSy  =  1  are  newly  the  same  as  those  found  in  the  more  resolved  case  of  Fig.  9 
(for  Dc  =  0.999,  and  =  0.016  cm).  For  each  flow  quantity,  the  contour  intervals  are 

chosen  to  be  the  same  in  Figs.  13a  and  13b.  The  effects  of  numerical  diffusion  on  the  flow  can 
then  be  inferred  by  examining  the  spreading  of  the  high-straln  region  at  the  interface  between 
the  two  streams  of  ^  based  on  the  distribution  and  density  of  contours,  and  by  examining  the 
numerical  stability  of  the  vortical  structures  based  on  the  vorticity  contours.  In  turn,  v*  -  based 
on  the  comparison  of  the  laminar  computed  and  theoretical  solutions  -  gives  us  an  indication 
of  the  amount  of  diffusion  responsible  for  that  spreading. 

In  the  case  of  the  calculation  on  grid  /  (aspect  ratio  unity),  numerical  diffusion  effects 
are  responsible  for  the  initial  spreading  of  the  step-velocity  profile  onto  an  inflectional  profile 
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confined  within  3-5  computational  cells,  and  are  otherwise  essentially  negligible  further  on 
downstream,  as  shown  in  Fig.  13a.  In  this  case,  the  numerical  diffusion  effects  are  isotropic, 
they  preserve  the  large  scale  features  of  the  flow  without  significant  distortion,  and  thus  mimic 
very  small  physical-viscooity  effects.  Although  the  general  features  of  the  flow  are  similar  in 
both  calculations,  it  is  apparent  that  Fig.  13b,  for  which  SxfSy  =  2.5,  shows  considerably 
larger  numerical  diffusion  effects.  This  is  indicated  by  the  significant  spreading  of  the  high- 
strain  regions  (particularly  on  the  faster  side),  and  by  the  increasingly  poor  resolution  of  the 
vortex  rolls  in  terms  of  vorticity  contours  as  we  move  downstream.  This  is  in  contrast  with 
the  results  in  Fig.  11  suggesting  that  no  significant  differences  should  be  expected  between 
calculations  on  grids  /  and  g.  major  difference,  however,  is  that  the  strain  direction  is  now  not 
necessarily  restricted  to  the  transverse  direction,  and  when  6x  >  6y  the  high-strain  region  is 
diffused  more  in  the  i-direction  than  in  the  y— direction. 
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4.  CONCLUSIONS 


We  have  shown  that  the  nonlinear  residual  numerical  diffusion  of  the  FCT  shear  flow  model 
with  1.  can  emulate  the  effects  of  physical  viscosity.  Effective  (global)  measures  of  the 

numerical  viscous  diffusion  in  the  model,  Ve,  were  obtained  by  comparing  the  numerical  spread 
of  low-Mach-number  simulated  mixing  layers  with  the  theoretical  spread  predicted  by  boundary 
layer  theory.  Associated  with  Vg,  the  profiles  of  calculated  z-momentum  at  different  streamwise 
locations  collapse  quite  well  in  the  vicinity  of  the  theoretical  results  when  plotted  as  a  function 
of  the  similarity  variable  rj.  We  found  that  the  global  numerical  diffusion  is  not  very  sensitive 
to  changes  in  free-stream  velocity  ratio  A,  and  can  be  reduced  to  a  desired  level  by  refining 
the  grid  spacing.  For  example,  for  Dg  =  0.999,  when  using  the  uniform  grid  a  in  Fig.  3,  and 
depending  on  A,  the  optimal  v*  was  found  to  vary  within  a  range  of  80%-100%  of  the  physical 
viscosity  value  for  air  at  the  temperature  and  pressure  conditions  of  the  simulations. 

When  non-square  uniform  grids  are  used,  the  numerical  diffusion  is  not  isotropic  and 
the  present  approach  for  its  global  interpretation  is  difficult  to  implement,  although  it  cam 
conceivably  be  used  to  obtain  bounds  on  its  magnitude.  We  found  that  the  large  scale  features 
of  the  shear  flow  can  be  significantly  distorted  by  the  effects  of  the  non-isotropic  numerical 
diffusion.  These  results  suggest  that  in  inviscid  simulations  of  shear  flows  where  numerical 
diffusion  is  the  only  cause  of  momentum  diffusion,  grids  with  aspect  ratios  as  close  to  unity  as 
possible  should  be  used  to  minimize  spurious  flow  distortions.  This  requirement  can  be  clearly 
relaxed  when  viscous  terms  are  included  in  the  model,  and  the  grid  size  is  chosen  to  ensure  that 
numerical  diffusion  is  negligible  as  compared  to  physical  diffusive  effects. 
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Fig.  1  a)  Flow  configuration  for  the  theoretical  study  of  the  laminar  spread  of  a  mixing  layer 
initialized  at  i  =  0  with  a  step-function  velocity  profile; 
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Fig.  2  Flow  configuration  for  the  simulated  mixing  layer. 
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Fig.  3  Computational  domains  and  grids  used  in  the  mixing  layer  simulations. 

•  ;  c  =  (U2  +  o)At/Ay  =  0.46. 

*  ;  20  evenly  spaced  computational  cells  are  used  in  the  y-direction  around  the  shear;  wall- 
boundaries  are  approached  through  geometrical  stretching  of  the  grid  (stretching  fac- 
tor=1.15). 

t  :  Gridding  in  x-direction  is  uniform. 

★  :  Grid  is  geometrically  stretched  in  i-dircction  for  i  >  4  (stretching  factor=1.03). 

o  :  c  =  0.4. 

1  :  60  evenly  spaced  cells  are  used  in  the  y-direction  around  the  shear;  wall-boundaries  are 
approached  through  geometrical  stretching  of  the  grid  (stretching  factor=1.15). 
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Fig.  4  Mixing  layer  simulation  on  grid  a  in  case  A  =  5,  Z?c  =  0.999 
deviations  and  <73. 
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Fig.  5  Mixing  layer  simulation  on  grid  a.  Dependence  of  diffusivities  on  each  side  of  the 
shear  layer  on  velocity  ratio  A  for  Dc  =  0.999.  The  arrows  indicate  least  values  of  Oi  at 
intersections  of  curves  of  partial  deviations  ai+  and  . 
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Fig.  8  Mixing  layer  simulation  on  grid  a  in  case  A  =  5,  De  =  0.990.  Comparison  of  momentum 
profiles  g{T})  for  u,  =  0.1 6cm* sec'*  at  selected  streamwise  locations  with  profile  {g{T})) 
predicted  by  boundary  layer  theory. 


Fig.  9  Mixing  layer  simulation  on  grids  o  -  c,  A  =  5.  Dependence  of  global  numerical  diffusion 
on  grid  spacing  for  =  0.990,  0.999.  The  dashed  lines  have  the  slope  Oi  =  1  predicted 
by  Eq.  (12). 
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Fig.  11  Mixing  layer  simulation  on  grids  a  and  d,  A  =  10,  =  0.999.  Dependence  of  globa 

numerical  diffusion  on  gridding  aspect  ratio. 


Fig.  12  Mixing  layer  simulation  on  grids  d  and  c,  A  =  10,  Dc  =  0.999.  Effects  of  grid  stretching 
in  x-direction  on  global  numerical  diffusion. 


Fig.  13  Fully  two-dimensional  mixing  layer  simulation  for  A  =  5.  Effects  of  gridding  aspect  ratio 
on  global  numerical  diffusion  for  Dc  =  0.999.  a)  6x  =  6y  (grid  /);  b)  6x  =  2.5^y  (grid 
g).  Instantaneous  flow  visualization  in  terms  of  contours  of  vorticity  (upper  frame)  and  of 
a  passive  scalar  convected  with  the  flow  velocity  (lower  frame).  The  contour  intervals  for 
each  flow  quantity  are  the  same  in  a)  and  b). 


